#ifndef ___DEFORM_H__
#define ___DEFORM_H__

#include <OGF/cells/common/common.h>
#include <OGF/cells/map/map.h>
#include <OGF/cells/map/map_attributes.h>

#include <OGF/basic/types/counted.h>
#include <OGF/math/attributes/attribute_adapter.h>
#include <OGF/surface/grob/surface.h>

//#include <OGF/math/linear_algebra/sparse_matrix.h>
#include <OGF/math/linear_algebra/vector.h>
#include <OGF/math/linear_algebra/matrix.h>
#include <OGF/math/geometry/types.h>


//#include <string>
#include <vector>
#include <set>
#include <iostream>

namespace OGF {

class Deform
{
// ================================== //
// ==== MISC ======================== //
// ================================== //
public:
	Deform();
	~Deform();

	void init(Surface *surface);		//Rebuild all the structures
	void lazy_init(Surface *surface);	//Init only if the surface has changed


private:
	Vector3d mult(const Matrix2d &m, const Vector3d &v) const;	//matrix-vector multiply
	Matrix2d fromAxisAngle(Vector3d axis, double angle);		//rotation matrix from axis & angle

	Surface *surface;

// ================================== //
// ==== REGION OPERATIONS =========== //
// ================================== //
public:

	/**
	* There are 3 important regions:
	*  BORDER: everything the user has frozen  (locked)
	*  GUARD: verticies adjacent to the border (laplace, and least-squares lock)
	*  CORE: the rest of the vertices in the ROI, except the anchor (laplace)
	*
	* Other sets:
	*  ANCHOR: just the anchor (locked, and laplace)
	*  ROI: union of BORDER, GUARD, CORE, and ANCHOR
	*  Everything Else: ignored outside the ROI
	*/

	//
	// Manage and query the border.
	///
	bool is_on_border(Map::Facet* face);
	void add_to_border(Map::Facet* face);
	void remove_from_border(Map::Facet* face);
	bool touches_border_face(Map::Vertex* vert);

	//
	// Manage and query the region of influence.
	// Used by the display code, and no one else?
	//
	bool is_in_region(Map::Vertex* vert);
	bool is_in_region(Map::Facet* facet);
	void grow_region();

	//
	// Manage and query the anchor vertex.
	//
	void set_anchor(Map::Vertex* vert);
	Map::Vertex* get_anchor();
	Vector3d get_true_anchor_normal();
	bool is_anchor_good();
private:
	void clear_region();
	void populate_verts();

	//TODO: support multiple anchor vertices?
	Map::Vertex* anchor;	// The anchor vertex (the one we drag around)
	
	enum CLASS{NIL=0, BORDER=1, GUARD=2, CORE=3, ANCHOR=4};
	MapFacetAttribute<bool>     *border;			//identify frozen faces
	MapVertexAttribute<CLASS>   *classify;			//classification (see enum CLASS)
	MapVertexAttribute<Point3d> *constraint_pos;	//original position of guard verticies

	// Three independent vertex sets.
	std::vector<Map::Vertex*> verts_border;	//locked
	std::vector<Map::Vertex*> verts_guard;	//laplace & soft-lock
	std::vector<Map::Vertex*> verts_core;	//laplace

	// Ordered union of above sets.
	std::vector< Map::Vertex* > verts;		// Verticies to affect. Order: [border, guard, core, anchor]
	MapVertexAttribute<int> *index;			// Reverse index into this list.

// ================================== //
// ==== DEFORMATION ================= //
// ================================== //
public:
	void bake_deformation();	//Saves the current mesh shape
	void do_fast_deformation();	//simple linear, for interaction
	void do_deformation();		//Full rotations and everything
	double gravity;

private:
	//
	// Gravity:
	//
	void calc_masses();
	void calc_areas();
	void calc_Laplace_weights();
	double calc_face_area(Map::Facet *f);
	void calc_alpha();
	void fill_tranpose_row(LinearSolver &solver, Map::Vertex *vert, int coord);
	void normalize_attribute(MapVertexAttribute<double> &att);
	
	//
	// Build the matrix system.
	//
	void compute_deltas(); //initial computation, to be done once on original mesh

	void do_laplace_solve();	//do the constrained laplace solve, after everything is setup.
	void build_system_solver(); //init the solver object. don't fill in the rows
	void setup_Laplace(int coord);
	void fill_constraint_row(Map::Vertex *vert, int coord);
	void fill_laplace_row(Map::Vertex *vert, int coord);

	bool has_laplace_row(Map::Vertex *v);
	bool has_constraint_row(Map::Vertex *v);

	//
	// Move coords between solver and surface, and vice-versa.
	//
	void solver_to_verts(int coord);
	void verts_to_solver(int coord);

	//
	// Rotation handling
	//
	void optimally_transform_deltas();
	void fill_neighbour_transform_rows(Map::Vertex *vert, LinearSolver *T_solver);

	//
	// Data
	//
	int GRAVITY_COORD;

	MapVertexAttribute<double> *mass;			// Mass of each vertex.
	std::vector<double> alpha;					// inv(L^T)*mass 

	MapFacetAttribute<double>  *area;			// Area of each face.
	MapVertexAttribute<double> *laplace_weight;	// Area-based weights for the laplace rows.

	MapVertexAttribute<Point3d>  *true_pos;		// Positions, prior to deformation
	MapVertexAttribute<Vector3d> *true_deltas;	// Delta vectors, prior to deformation
	MapVertexAttribute<Vector3d> *deltas;		// rotated Laplacian delta vectors - use these ones in solver

	LinearSolver *solver;						// The linear solver

	//
	// Paramaters:
	//
	double constraint_weight;
};

extern Deform deform;

} //namespace


#endif // ___DEFORM_H__
